### REPLICATION FILE -- ANALYSIS
### Jonathan Homola
### "Partisanship and Perceived Threats about Immigration"
### Party Politics

rm(list=ls())
set.seed(12435)
options(stringsAsFactors=F)
library(stargazer); library(xtable); library(psych); library(AER)

data <- read.csv("immig_analysis.csv")

data$treat[is.na(data$aug2015wt1)] <- NA               ## T: 788, C: 750, N: 1671
data1 <- data[(data$manip2!=0 | is.na(data$manip2)),]  ## T: 732, C: 750, N: 1615
data2 <- data[(data$manip1!=0 | is.na(data$manip1)),]  ## T: 439, C: 750, N: 1322
names(data)

## Start calculating t-tests of treatment effects

## Overall
result1 <- pvals1 <- NULL
confint1 <- matrix(NA, ncol=2, nrow=6)
for (i in 6:11){
  temp <- t.test(as.numeric(data[,i]) ~ data$treat)
  result1 <- c(result1, temp$estimate[[2]], temp$estimate[[1]], (temp$estimate[[2]] - temp$estimate[[1]]))
  pvals1 <- c(pvals1, temp$p.value)
  confint1[12-i,] <- -temp$conf.int[c(1,2)]
}
result1 <- c(result1, nrow(data[!is.na(data$order_viol_comm),]))

## Democrats only
dems <- subset(data, PID3==1)
for (i in 6:11){
  temp <- t.test(as.numeric(dems[,i]) ~ dems$treat)
  result1 <- c(result1, temp$estimate[[2]], temp$estimate[[1]], (temp$estimate[[2]] - temp$estimate[[1]]))
  pvals1 <- c(pvals1, temp$p.value)
}
result1 <- c(result1, nrow(dems[!is.na(dems$order_viol_comm),]))

## Republicans only
reps <- subset(data, PID3==2)
for (i in 6:11){
  temp <- t.test(as.numeric(reps[,i]) ~ reps$treat)
  result1 <- c(result1, temp$estimate[[2]], temp$estimate[[1]], (temp$estimate[[2]] - temp$estimate[[1]]))
  pvals1 <- c(pvals1, temp$p.value)
}
result1 <- c(result1, nrow(reps[!is.na(reps$order_viol_comm),]))

## Rep-Dem diff
treat <- subset(data, (treat==1 & PID3<3))
treat$dem <- ifelse(treat$PID3==1, 1, 0)
cont <- subset(data, (treat==0 & PID3<3))
cont$dem <- ifelse(cont$PID3==1, 1, 0)
data_demrep <- subset(data, PID3<3)
data_demrep$rep <- ifelse(data_demrep$PID3==2, 1, 0)
for (i in 6:11){
  temp <- t.test(as.numeric(treat[,i]) ~ treat$dem)
  result1 <- c(result1, temp$estimate[[1]] - temp$estimate[[2]])
  pvals1 <- c(pvals1, temp$p.value)
  temp <- t.test(as.numeric(cont[,i]) ~ cont$dem)
  result1 <- c(result1, temp$estimate[[1]] - temp$estimate[[2]])
  pvals1 <- c(pvals1, temp$p.value)
  result1 <- c(result1, summary(lm(data_demrep[,i] ~ data_demrep$treat*data_demrep$rep))$coef[4,1])
  pvals1 <- c(pvals1, summary(lm(data_demrep[,i] ~ data_demrep$treat*data_demrep$rep))$coef[4,4])
}


## Create results/output tables (Table 2)
results <- matrix(round(result1, 2), ncol=4, nrow=19, byrow=F)
results
xtable(results)
results2 <- data.frame(cbind(c(rep(c("Treatment", "Control", "ATE"), 6), "N"), results))
print(xtable(results2), include.rownames = FALSE)

## Check p-values
round(pvals1, 3)[1:6]
round(pvals1, 3)[7:12]
round(pvals1, 3)[13:18]
round(pvals1, 3)[19:36]

## create plot 1 (overall treatment effects)
catNames <- c("American culture", "National identity",
              "Economy, national", "Economy, household",
              "Violence, national", "Violence, community")
#pdf("treat1.pdf", height=6,width=8)
par(oma=c(1,3.5,0.3,1), mar=c(3,3.5,0.3,1))
plot(0,0, type="n", xaxt="n", yaxt="n", xlab="", ylab="", main="",
     xlim=c(-0.15, 0.25), ylim=c(0,6))
abline(h=1:5, col="gray60", lty=2)

# plot treatment effects, confidence intervals and null effect
points(result1[c(18,15,12,9,6,3)], seq(0.5,5.5,1), pch=16)
segments(confint1[,1], seq(0.5,5.5,1), confint1[,2], seq(0.5,5.5,1))
abline(v=0, col="red")

# plot axes
axis(1, tck=0, cex.axis=0.8, cex.lab=0.8, mgp=c(0.3, 0.3, 0), lty=0)
axis(2, tck=0, cex.axis=0.8, cex.lab=0.8, at = seq(0.5,5.5,1), labels=catNames,
     mgp=c(0.3, 0.3, 0), lty=0, las=1)
title(xlab = "Treatment Effect", line = 1.5, cex.lab = 1, ylab="")
#dev.off()


#####
##### APPENDIX
#####

## SM4: Randomization/balance checks
## make Independents baseline category for PID
table(data$PID3) ## 1=D, 2=R, 3=I
table(data$ppreg4) ## 1=NE, 2=MW, 3=S, 4=W
data$PID3 <- relevel(as.factor(data$PID3), ref=3)
data$ppeducat <- factor(data$ppeducat, levels = c("Less than high school", "High school", "Some college", "Bachelor's degree or higher"))
table(data$ppreg4) ## 1=NE, 2=MW, 3=S, 4=W
rand1 <- glm(treat ~ as.numeric(female), data, family="binomial")
rand2 <- glm(treat ~ income6, data, family="binomial")
rand3 <- glm(treat ~ as.factor(PID3), data, family="binomial")
rand4 <- glm(treat ~ as.factor(ppeducat), data, family="binomial")
rand5 <- glm(treat ~ as.factor(ppreg4), data, family="binomial")
rand6 <- glm(treat ~ as.numeric(female) + income6 + as.factor(PID3) +
               as.factor(ppeducat) + as.factor(ppreg4), data, family="binomial")

# Table SM4.1
stargazer(rand1, rand2, rand3, rand4, rand5, rand6,
          omit.stat=c("f", "ser", "adj.rsq"), star.cutoffs = c(0.05, NA, NA),
          column.sep.width="0.5pt", digits=3, no.space=T,
          dep.var.caption="Outcome Variable: Treatment Assignment", 
          dep.var.labels.include = FALSE)

## Randomization/balance checks: DEMOCRATS
dems <- subset(data, PID3==1)
rand7 <- glm(treat ~ as.numeric(female), dems, family="binomial")
rand8 <- glm(treat ~ as.numeric(income6), dems, family="binomial")
rand9 <- glm(treat ~ as.factor(ppeducat), dems, family="binomial")
rand10 <- glm(treat ~ as.factor(ppreg4), dems, family="binomial")
rand11 <- glm(treat ~ as.numeric(female) + as.numeric(income6) +
                as.factor(ppeducat) + as.factor(ppreg4), dems, family="binomial")

# Table SM4.2
stargazer(rand7, rand8, rand9, rand10, rand11,
          omit.stat=c("f", "ser", "adj.rsq"), star.cutoffs = c(0.05, NA, NA),
          column.sep.width="0.5pt", digits=3, no.space=T,
          dep.var.caption="Outcome Variable: Treatment Assignment", 
          dep.var.labels.include = FALSE)

## Randomization/balance checks: REPUBLICANS
reps <- subset(data, PID3==2)
rand12 <- glm(treat ~ as.numeric(female), reps, family="binomial")
rand13 <- glm(treat ~ as.numeric(income6), reps, family="binomial")
rand14 <- glm(treat ~ as.factor(ppeducat), reps, family="binomial")
rand15 <- glm(treat ~ as.factor(ppreg4), reps, family="binomial")
rand16 <- glm(treat ~ as.numeric(female) + as.numeric(income6) +
                as.factor(ppeducat) + as.factor(ppreg4), reps, family="binomial")

# Table SM4.3
stargazer(rand12, rand13, rand14, rand15, rand16,
          omit.stat=c("f", "ser", "adj.rsq"), star.cutoffs = c(0.05, NA, NA),
          column.sep.width="0.5pt", digits=3, no.space=T,
          dep.var.caption="Outcome Variable: Treatment Assignment", 
          dep.var.labels.include = FALSE)


## SM3: Descriptive statistics
data$female <- as.numeric(data$female)-1
data$HighSchool <- ifelse(as.numeric(data$ppeducat)==2, 1, 0)
data$SomeCollege <- ifelse(as.numeric(data$ppeducat)==3, 1, 0)
data$BA <- ifelse(as.numeric(data$ppeducat)==4, 1, 0)
data$Midwest <- ifelse(as.numeric(data$ppreg4)==2, 1, 0)
data$South <- ifelse(as.numeric(data$ppreg4)==3, 1, 0)
data$West <- ifelse(as.numeric(data$ppreg4)==4, 1, 0)
data$Democrat <- ifelse(data$PID3==1, 1, 0)
data$Democrat <- ifelse(is.na(data$PID3), 0, data$Democrat)
data$Republican <- ifelse(data$PID3==2, 1, 0)
data$Republican <- ifelse(is.na(data$PID3), 0, data$Republican)
data$White <- ifelse(data$ppethm==1, 1, 0)

desc1 <- data[!is.na(data$treat),]
desc2 <- subset(data, treat==1)
desc3 <- subset(data, treat==0)

# Table SM3.1
xtable(describe(desc1[,c("female", "AGE_2016", "income6", "White", "HighSchool", "SomeCollege", "BA",
                         "Midwest", "South", "West", "Democrat", "Republican")])[,c(2,3,4,8,9)], digits=3)

xtable(describe(desc2[,c("female", "AGE_2016", "income6", "White", "HighSchool", "SomeCollege", "BA",
                         "Midwest", "South", "West", "Democrat", "Republican")])[,c(2,3,4,8,9)], digits=3)

xtable(describe(desc3[,c("female", "AGE_2016", "income6", "White", "HighSchool", "SomeCollege", "BA",
                         "Midwest", "South", "West", "Democrat", "Republican")])[,c(2,3,4,8,9)], digits=3)

# Table SM3.2
xtable(describe(desc1[,c("viol_comm", "viol_nat", "econ_hh", "econ_nat", "ident",
                         "culture")])[,c(2,3,4,8,9)], digits=3)

xtable(describe(desc2[,c("viol_comm", "viol_nat", "econ_hh", "econ_nat", "ident",
                         "culture")])[,c(2,3,4,8,9)], digits=3)

xtable(describe(desc3[,c("viol_comm", "viol_nat", "econ_hh", "econ_nat", "ident",
                         "culture")])[,c(2,3,4,8,9)], digits=3)


## SM6: Robustness: Income/Education table
## Start calculating t-tests of treatment effects
## Overall
result1 <- pvals1 <- NULL
for (i in 6:11){
  temp <- t.test(as.numeric(data[,i]) ~ data$treat)
  result1 <- c(result1, temp$estimate[[2]], temp$estimate[[1]], (temp$estimate[[2]] - temp$estimate[[1]]))
  pvals1 <- c(pvals1, temp$p.value)
}
result1 <- c(result1, nrow(data[!is.na(data$order_viol_comm),]))

## Income below $50,000
lowinc <- subset(data, as.numeric(income6)<4)
for (i in 6:11){
  temp <- t.test(as.numeric(lowinc[,i]) ~ lowinc$treat)
  result1 <- c(result1, temp$estimate[[2]], temp$estimate[[1]], (temp$estimate[[2]] - temp$estimate[[1]]))
  pvals1 <- c(pvals1, temp$p.value)
}
result1 <- c(result1, nrow(lowinc[!is.na(lowinc$order_viol_comm),]))

## Income above $50,000
hiinc <- subset(data, as.numeric(income6)>3)
for (i in 6:11){
  temp <- t.test(as.numeric(hiinc[,i]) ~ hiinc$treat)
  result1 <- c(result1, temp$estimate[[2]], temp$estimate[[1]], (temp$estimate[[2]] - temp$estimate[[1]]))
  pvals1 <- c(pvals1, temp$p.value)
}
result1 <- c(result1, nrow(hiinc[!is.na(hiinc$order_viol_comm),]))

## Education less than HS, HS, some college
lowedu <- subset(data, as.numeric(ppeducat)<4)
for (i in 6:11){
  temp <- t.test(as.numeric(lowedu[,i]) ~ lowedu$treat)
  result1 <- c(result1, temp$estimate[[2]], temp$estimate[[1]], (temp$estimate[[2]] - temp$estimate[[1]]))
  pvals1 <- c(pvals1, temp$p.value)
}
result1 <- c(result1, nrow(lowedu[!is.na(lowedu$order_viol_comm),]))

## Education Bachelor's degree or higher
hiedu <- subset(data, as.numeric(ppeducat)>3)
for (i in 6:11){
  temp <- t.test(as.numeric(hiedu[,i]) ~ hiedu$treat)
  result1 <- c(result1, temp$estimate[[2]], temp$estimate[[1]], (temp$estimate[[2]] - temp$estimate[[1]]))
  pvals1 <- c(pvals1, temp$p.value)
}
result1 <- c(result1, nrow(hiedu[!is.na(hiedu$order_viol_comm),]))

## Create results/output tables
results <- matrix(round(result1, 2), ncol=5, nrow=19, byrow=F)
results
xtable(results)
results2 <- data.frame(cbind(c(rep(c("Treatment", "Control", "Difference"), 6), "N"), results))

# Table SM6.1
print(xtable(results2), include.rownames = FALSE)

## Check p-values
round(pvals1, 3)[1:6]
round(pvals1, 3)[7:12]
round(pvals1, 3)[13:18]
round(pvals1, 3)[19:24]
round(pvals1, 3)[25:30]

## ROBUSTNESS: WHITES ONLY
data$White <- ifelse(data$ppethm==1, 1, 0)

## Start calculating t-tests of treatment effects
datawh <- subset(data, White==1)
datanw <- subset(data, White==0)

## Overall -- whites only
result1 <- pvals1 <- NULL
confint1 <- matrix(NA, ncol=2, nrow=6)
for (i in 6:11){
  temp <- t.test(as.numeric(datawh[,i]) ~ datawh$treat)
  result1 <- c(result1, temp$estimate[[2]], temp$estimate[[1]], (temp$estimate[[2]] - temp$estimate[[1]]))
  pvals1 <- c(pvals1, temp$p.value)
  confint1[12-i,] <- -temp$conf.int[c(1,2)]
}
result1 <- c(result1, nrow(datawh[!is.na(datawh$order_viol_comm),]))

## Nonwhites only
for (i in 6:11){
  temp <- t.test(as.numeric(datanw[,i]) ~ datanw$treat)
  result1 <- c(result1, temp$estimate[[2]], temp$estimate[[1]], (temp$estimate[[2]] - temp$estimate[[1]]))
  pvals1 <- c(pvals1, temp$p.value)
}
result1 <- c(result1, nrow(datanw[!is.na(datanw$order_viol_comm),]))

## White Democrats only
dems <- subset(datawh, PID3==1)
for (i in 6:11){
  temp <- t.test(as.numeric(dems[,i]) ~ dems$treat)
  result1 <- c(result1, temp$estimate[[2]], temp$estimate[[1]], (temp$estimate[[2]] - temp$estimate[[1]]))
  pvals1 <- c(pvals1, temp$p.value)
}
result1 <- c(result1, nrow(dems[!is.na(dems$order_viol_comm),]))

## White Republicans only
reps <- subset(datawh, PID3==2)
for (i in 6:11){
  temp <- t.test(as.numeric(reps[,i]) ~ reps$treat)
  result1 <- c(result1, temp$estimate[[2]], temp$estimate[[1]], (temp$estimate[[2]] - temp$estimate[[1]]))
  pvals1 <- c(pvals1, temp$p.value)
}
result1 <- c(result1, nrow(reps[!is.na(reps$order_viol_comm),]))

## Create results/output tables
results <- matrix(round(result1, 2), ncol=4, nrow=19, byrow=F)
results
xtable(results)
results2 <- data.frame(cbind(c(rep(c("Treatment", "Control", "Difference"), 6), "N"), results))

# Table SM6.2
print(xtable(results2), include.rownames = FALSE)

## Check p-values
round(pvals1, 3)[1:6]
round(pvals1, 3)[7:12]
round(pvals1, 3)[13:18]
round(pvals1, 3)[19:24]


## Figure SM6.1
## create plot  (overall treatment effects)
catNames <- c("American culture", "National identity",
              "Economy, national", "Economy, household",
              "Violence, national", "Violence, community")
#pdf("treat1_white.pdf", height=6,width=8)
par(oma=c(1,3.5,0.3,1), mar=c(3,3.5,0.3,1))
plot(0,0, type="n", xaxt="n", yaxt="n", xlab="", ylab="", main="",
     xlim=c(-0.15, 0.25), ylim=c(0,6))
abline(h=1:5, col="gray60", lty=2)

# plot treatment effects, confidence intervals and null effect
points(result1[c(18,15,12,9,6,3)], seq(0.5,5.5,1), pch=16)
segments(confint1[,1], seq(0.5,5.5,1), confint1[,2], seq(0.5,5.5,1))
abline(v=0, col="red")

# plot axes
axis(1, tck=0, cex.axis=0.8, cex.lab=0.8, mgp=c(0.3, 0.3, 0), lty=0)
axis(2, tck=0, cex.axis=0.8, cex.lab=0.8, at = seq(0.5,5.5,1), labels=catNames,
     mgp=c(0.3, 0.3, 0), lty=0, las=1)
title(xlab = "Treatment Effect", line = 1.5, cex.lab = 1, ylab="")
#dev.off()


#### ROBUSTNESS: HIGH/LOW IMMIGRATION STATES
census <- read.csv("ACS_13_1YR_CP02_with_ann.csv")
states <- census$GEO.display.label
foreign <- census$HC01_VC136
census <- data.frame(state=states[2:52], foreign=foreign[2:52])
census$foreign <- as.numeric(census$foreign)
census$st <- state.abb[match(census$state,state.name)]
census$st[census$state=="District of Columbia"] <- "DC"
summary(census$foreign)   # Median share of foreign-born pop: 6.8
census$highshare <- ifelse(census$foreign>6.8, 1, 0)
census1 <- data.frame(ppstate=census$st, share=census$foreign, highshare=census$highshare)
data1 <- merge(data, census1, by="ppstate")

mod1 <- lm(viol_comm ~ treat, data1)
mod2 <- lm(viol_nat ~ treat, data1)
mod3 <- lm(econ_hh ~ treat, data1)
mod4 <- lm(econ_nat ~ treat, data1)
mod5 <- lm(ident ~ treat, data1)
mod6 <- lm(culture ~ treat, data1)

mod1a <- lm(viol_comm ~ treat + share, data1)
mod2a <- lm(viol_nat ~ treat + share, data1)
mod3a <- lm(econ_hh ~ treat + share, data1)
mod4a <- lm(econ_nat ~ treat + share, data1)
mod5a <- lm(ident ~ treat + share, data1)
mod6a <- lm(culture ~ treat + share, data1)

# Table SM6.3 -- first panel (OVERALL)
stargazer(mod1, mod2, mod3, mod4, mod5, mod6)
stargazer(mod1a, mod2a, mod3a, mod4a, mod5a, mod6a)

dems <- subset(data1, PID3==1)
mod1 <- lm(viol_comm ~ treat, dems)
mod2 <- lm(viol_nat ~ treat, dems)
mod3 <- lm(econ_hh ~ treat, dems)
mod4 <- lm(econ_nat ~ treat, dems)
mod5 <- lm(ident ~ treat, dems)
mod6 <- lm(culture ~ treat, dems)

mod1a <- lm(viol_comm ~ treat + share, dems)
mod2a <- lm(viol_nat ~ treat + share, dems)
mod3a <- lm(econ_hh ~ treat + share, dems)
mod4a <- lm(econ_nat ~ treat + share, dems)
mod5a <- lm(ident ~ treat + share, dems)
mod6a <- lm(culture ~ treat + share, dems)

# Table SM6.3 -- second panel (DEMOCRATS)
stargazer(mod1, mod2, mod3, mod4, mod5, mod6)
stargazer(mod1a, mod2a, mod3a, mod4a, mod5a, mod6a)

reps <- subset(data1, PID3==2)
mod1 <- lm(viol_comm ~ treat, reps)
mod2 <- lm(viol_nat ~ treat, reps)
mod3 <- lm(econ_hh ~ treat, reps)
mod4 <- lm(econ_nat ~ treat, reps)
mod5 <- lm(ident ~ treat, reps)
mod6 <- lm(culture ~ treat, reps)

mod1a <- lm(viol_comm ~ treat + share, reps)
mod2a <- lm(viol_nat ~ treat + share, reps)
mod3a <- lm(econ_hh ~ treat + share, reps)
mod4a <- lm(econ_nat ~ treat + share, reps)
mod5a <- lm(ident ~ treat + share, reps)
mod6a <- lm(culture ~ treat + share, reps)

# Table SM6.3 -- third panel (REPUBLICANS)
stargazer(mod1, mod2, mod3, mod4, mod5, mod6)
stargazer(mod1a, mod2a, mod3a, mod4a, mod5a, mod6a)


## SM5: Attention checks
### LATE
table(data$treat, useNA = "ifany")
table(data$manip1, useNA = "ifany")
table(data$manip2, useNA = "ifany")
table(data$treat, data$manip1, useNA = "ifany")
data$manip1[data$treat==0] <- 0
data$manip2[data$treat==0] <- 0

dems <- subset(data, PID3==1)
reps <- subset(data, PID3==2)

## MONTANA OR MISSISSIPPI
iv1 <- ivreg(viol_comm~manip2 | treat, data=data)
iv2 <- ivreg(viol_nat~manip2 | treat, data=data)
iv3 <- ivreg(econ_hh~manip2 | treat, data=data)
iv4 <- ivreg(econ_nat~manip2 | treat, data=data)
iv5 <- ivreg(ident~manip2 | treat, data=data)
iv6 <- ivreg(culture~manip2 | treat, data=data)

iv7 <- ivreg(viol_comm~manip2 | treat, data=dems)
iv8 <- ivreg(viol_nat~manip2 | treat, data=dems)
iv9 <- ivreg(econ_hh~manip2 | treat, data=dems)
iv10 <- ivreg(econ_nat~manip2 | treat, data=dems)
iv11 <- ivreg(ident~manip2 | treat, data=dems)
iv12 <- ivreg(culture~manip2 | treat, data=dems)

iv13 <- ivreg(viol_comm~manip2 | treat, data=reps)
iv14 <- ivreg(viol_nat~manip2 | treat, data=reps)
iv15 <- ivreg(econ_hh~manip2 | treat, data=reps)
iv16 <- ivreg(econ_nat~manip2 | treat, data=reps)
iv17 <- ivreg(ident~manip2 | treat, data=reps)
iv18 <- ivreg(culture~manip2 | treat, data=reps)

result1 <- c(sum(iv1$coefficients), iv1$coefficients[[1]], iv1$coefficients[[2]],
             sum(iv2$coefficients), iv2$coefficients[[1]], iv2$coefficients[[2]],
             sum(iv3$coefficients), iv3$coefficients[[1]], iv3$coefficients[[2]],
             sum(iv4$coefficients), iv4$coefficients[[1]], iv4$coefficients[[2]],
             sum(iv5$coefficients), iv5$coefficients[[1]], iv5$coefficients[[2]],
             sum(iv6$coefficients), iv6$coefficients[[1]], iv6$coefficients[[2]], 1538,
             sum(iv7$coefficients), iv7$coefficients[[1]], iv7$coefficients[[2]],
             sum(iv8$coefficients), iv8$coefficients[[1]], iv8$coefficients[[2]],
             sum(iv9$coefficients), iv9$coefficients[[1]], iv9$coefficients[[2]],
             sum(iv10$coefficients), iv10$coefficients[[1]], iv10$coefficients[[2]], 
             sum(iv11$coefficients), iv11$coefficients[[1]], iv11$coefficients[[2]], 
             sum(iv12$coefficients), iv12$coefficients[[1]], iv12$coefficients[[2]], 767,
             sum(iv13$coefficients), iv13$coefficients[[1]], iv13$coefficients[[2]], 
             sum(iv14$coefficients), iv14$coefficients[[1]], iv14$coefficients[[2]], 
             sum(iv15$coefficients), iv15$coefficients[[1]], iv15$coefficients[[2]],
             sum(iv16$coefficients), iv16$coefficients[[1]], iv16$coefficients[[2]], 
             sum(iv17$coefficients), iv17$coefficients[[1]], iv17$coefficients[[2]],
             sum(iv18$coefficients), iv18$coefficients[[1]], iv18$coefficients[[2]], 651)

pvals <- c(summary(iv1)$coefficients[2,4], summary(iv2)$coefficients[2,4],
           summary(iv3)$coefficients[2,4], summary(iv4)$coefficients[2,4],
           summary(iv5)$coefficients[2,4], summary(iv6)$coefficients[2,4],
           summary(iv7)$coefficients[2,4], summary(iv8)$coefficients[2,4],
           summary(iv9)$coefficients[2,4], summary(iv10)$coefficients[2,4],
           summary(iv11)$coefficients[2,4], summary(iv12)$coefficients[2,4],
           summary(iv13)$coefficients[2,4], summary(iv14)$coefficients[2,4],
           summary(iv15)$coefficients[2,4], summary(iv16)$coefficients[2,4],
           summary(iv17)$coefficients[2,4], summary(iv18)$coefficients[2,4])

## Create results/output tables
results <- matrix(round(result1, 2), ncol=3, nrow=19, byrow=F)
results
xtable(results)
results2 <- data.frame(cbind(c(rep(c("Treatment", "Control", "Difference"), 6), "N"), results))

# Table SM5.1
print(xtable(results2), include.rownames = FALSE)

## Check p-values
round(pvals, 3)[1:6]
round(pvals, 3)[7:12]
round(pvals, 4)[13:18]

## ONLY MONTANA
iv1 <- ivreg(viol_comm~manip1 | treat, data=data)
iv2 <- ivreg(viol_nat~manip1 | treat, data=data)
iv3 <- ivreg(econ_hh~manip1 | treat, data=data)
iv4 <- ivreg(econ_nat~manip1 | treat, data=data)
iv5 <- ivreg(ident~manip1 | treat, data=data)
iv6 <- ivreg(culture~manip1 | treat, data=data)

iv7 <- ivreg(viol_comm~manip1 | treat, data=dems)
iv8 <- ivreg(viol_nat~manip1 | treat, data=dems)
iv9 <- ivreg(econ_hh~manip1 | treat, data=dems)
iv10 <- ivreg(econ_nat~manip1 | treat, data=dems)
iv11 <- ivreg(ident~manip1 | treat, data=dems)
iv12 <- ivreg(culture~manip1 | treat, data=dems)

iv13 <- ivreg(viol_comm~manip1 | treat, data=reps)
iv14 <- ivreg(viol_nat~manip1 | treat, data=reps)
iv15 <- ivreg(econ_hh~manip1 | treat, data=reps)
iv16 <- ivreg(econ_nat~manip1 | treat, data=reps)
iv17 <- ivreg(ident~manip1 | treat, data=reps)
iv18 <- ivreg(culture~manip1 | treat, data=reps)

result1 <- c(sum(iv1$coefficients), iv1$coefficients[[1]], iv1$coefficients[[2]],
             sum(iv2$coefficients), iv2$coefficients[[1]], iv2$coefficients[[2]],
             sum(iv3$coefficients), iv3$coefficients[[1]], iv3$coefficients[[2]],
             sum(iv4$coefficients), iv4$coefficients[[1]], iv4$coefficients[[2]],
             sum(iv5$coefficients), iv5$coefficients[[1]], iv5$coefficients[[2]],
             sum(iv6$coefficients), iv6$coefficients[[1]], iv6$coefficients[[2]], 1538,
             sum(iv7$coefficients), iv7$coefficients[[1]], iv7$coefficients[[2]],
             sum(iv8$coefficients), iv8$coefficients[[1]], iv8$coefficients[[2]],
             sum(iv9$coefficients), iv9$coefficients[[1]], iv9$coefficients[[2]],
             sum(iv10$coefficients), iv10$coefficients[[1]], iv10$coefficients[[2]], 
             sum(iv11$coefficients), iv11$coefficients[[1]], iv11$coefficients[[2]], 
             sum(iv12$coefficients), iv12$coefficients[[1]], iv12$coefficients[[2]], 767,
             sum(iv13$coefficients), iv13$coefficients[[1]], iv13$coefficients[[2]], 
             sum(iv14$coefficients), iv14$coefficients[[1]], iv14$coefficients[[2]], 
             sum(iv15$coefficients), iv15$coefficients[[1]], iv15$coefficients[[2]],
             sum(iv16$coefficients), iv16$coefficients[[1]], iv16$coefficients[[2]], 
             sum(iv17$coefficients), iv17$coefficients[[1]], iv17$coefficients[[2]],
             sum(iv18$coefficients), iv18$coefficients[[1]], iv18$coefficients[[2]], 651)

pvals <- c(summary(iv1)$coefficients[2,4], summary(iv2)$coefficients[2,4],
           summary(iv3)$coefficients[2,4], summary(iv4)$coefficients[2,4],
           summary(iv5)$coefficients[2,4], summary(iv6)$coefficients[2,4],
           summary(iv7)$coefficients[2,4], summary(iv8)$coefficients[2,4],
           summary(iv9)$coefficients[2,4], summary(iv10)$coefficients[2,4],
           summary(iv11)$coefficients[2,4], summary(iv12)$coefficients[2,4],
           summary(iv13)$coefficients[2,4], summary(iv14)$coefficients[2,4],
           summary(iv15)$coefficients[2,4], summary(iv16)$coefficients[2,4],
           summary(iv17)$coefficients[2,4], summary(iv18)$coefficients[2,4])

## Create results/output tables
results <- matrix(round(result1, 2), ncol=3, nrow=19, byrow=F)
results
xtable(results)
results2 <- data.frame(cbind(c(rep(c("Treatment", "Control", "Difference"), 6), "N"), results))

# Table SM5.2
print(xtable(results2), include.rownames = FALSE)

## Check p-values
round(pvals, 3)[1:6]
round(pvals, 3)[7:12]
round(pvals, 4)[13:18]


## SM6: Placebo information treatment
mturk <- read.csv("mturk_immig.csv")

table(mturk$place, useNA="ifany") ## 439 control, 379 placebo

mturk$dem <- ifelse(mturk$Q5=="Democrat", 1, NA)
mturk$dem <- ifelse(mturk$Q6=="Democratic Party", 1, mturk$dem)
mturk$dem <- ifelse(mturk$Q5=="Republican", 0, mturk$dem)
mturk$dem <- ifelse(mturk$Q6=="Republican Party", 0, mturk$dem)

## 1 vio comm, 2 vio count
## 3 econ hh, 4 econ count
## 5 identity, 6 culture
mturk$vio1 <- ifelse(mturk$Q16_1=="Agree strongly", 4, NA)
mturk$vio1 <- ifelse(mturk$Q16_1=="Agree somewhat", 3, mturk$vio1)
mturk$vio1 <- ifelse(mturk$Q16_1=="Disagree somewhat", 2, mturk$vio1)
mturk$vio1 <- ifelse(mturk$Q16_1=="Disagree strongly", 1, mturk$vio1)

mturk$vio2 <- ifelse(mturk$Q16_2=="Agree strongly", 4, NA)
mturk$vio2 <- ifelse(mturk$Q16_2=="Agree somewhat", 3, mturk$vio2)
mturk$vio2 <- ifelse(mturk$Q16_2=="Disagree somewhat", 2, mturk$vio2)
mturk$vio2 <- ifelse(mturk$Q16_2=="Disagree strongly", 1, mturk$vio2)

mturk$eco1 <- ifelse(mturk$Q16_3=="Agree strongly", 4, NA)
mturk$eco1 <- ifelse(mturk$Q16_3=="Agree somewhat", 3, mturk$eco1)
mturk$eco1 <- ifelse(mturk$Q16_3=="Disagree somewhat", 2, mturk$eco1)
mturk$eco1 <- ifelse(mturk$Q16_3=="Disagree strongly", 1, mturk$eco1)

mturk$eco2 <- ifelse(mturk$Q16_4=="Agree strongly", 4, NA)
mturk$eco2 <- ifelse(mturk$Q16_4=="Agree somewhat", 3, mturk$eco2)
mturk$eco2 <- ifelse(mturk$Q16_4=="Disagree somewhat", 2, mturk$eco2)
mturk$eco2 <- ifelse(mturk$Q16_4=="Disagree strongly", 1, mturk$eco2)

mturk$cul1 <- ifelse(mturk$Q16_5=="Agree strongly", 4, NA)
mturk$cul1 <- ifelse(mturk$Q16_5=="Agree somewhat", 3, mturk$cul1)
mturk$cul1 <- ifelse(mturk$Q16_5=="Disagree somewhat", 2, mturk$cul1)
mturk$cul1 <- ifelse(mturk$Q16_5=="Disagree strongly", 1, mturk$cul1)

mturk$cul2 <- ifelse(mturk$Q16_6=="Agree strongly", 4, NA)
mturk$cul2 <- ifelse(mturk$Q16_6=="Agree somewhat", 3, mturk$cul2)
mturk$cul2 <- ifelse(mturk$Q16_6=="Disagree somewhat", 2, mturk$cul2)
mturk$cul2 <- ifelse(mturk$Q16_6=="Disagree strongly", 1, mturk$cul2)

dems <- subset(mturk, dem==1)
reps <- subset(mturk, dem==0)
table(dems$place, useNA="ifany") ## 257 control, 238 placebo
table(reps$place, useNA="ifany") ## 180 control, 141 placebo

### create plot for placebo vs control
mturk_placebo <- mturk
result1 <- pvals1 <- NULL
confint1 <- matrix(NA, ncol=2, nrow=6)
for (i in 50:55){
  temp <- t.test(as.numeric(mturk_placebo[,i]) ~ mturk_placebo$place)
  result1 <- c(result1, temp$estimate[[2]], temp$estimate[[1]], ((temp$estimate[[2]]) - (temp$estimate[[1]])))
  pvals1 <- c(pvals1, temp$p.value)
  confint1[56-i,] <- (-1)*temp$conf.int[c(1,2)]
}

## Figure SM6.2
## create plot
catNames <- c("American culture", "National identity",
              "Economy, national", "Economy, household",
              "Violence, national", "Violence, community")
#pdf("placebo1.pdf", height=6,width=8)
par(oma=c(1,3.5,0.3,1), mar=c(3,3.5,0.3,1))
plot(0,0, type="n", xaxt="n", yaxt="n", xlab="", ylab="", main="",
     xlim=c(-0.35, 0.25), ylim=c(0,6))
abline(h=1:5, col="gray60", lty=2)

# plot treatment effects, confidence intervals and null effect
points(result1[c(18,15,12,9,6,3)], seq(0.5,5.5,1), pch=16)
segments(confint1[,1], seq(0.5,5.5,1), confint1[,2], seq(0.5,5.5,1))
abline(v=0, col="red")

# plot axes
axis(1, tck=0, cex.axis=0.8, cex.lab=0.8, mgp=c(0.3, 0.3, 0), lty=0)
axis(2, tck=0, cex.axis=0.8, cex.lab=0.8, at = seq(0.5,5.5,1), labels=catNames,
     mgp=c(0.3, 0.3, 0), lty=0, las=1)
title(xlab = "Placebo Effect", line = 1.5, cex.lab = 1, ylab="")
#dev.off()